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Abstract. New interlayer intermolecular potential model was proposed and it rep- 
resented "ABAB" staking of graphite. Hydrogen atom sputtering on graphite surface 
was investigated using molecular dynamics simulation. In the initial short time period, 
maintaining the flat structure of graphenes, hydrogen atoms brought about the differ- 
ence interaction process in each incident energy. The first graphene often adsorbed 
5 eV hydrogen atoms and reflected almost all of 15 eV hydrogen atoms. The hydrogen 
atoms which were injected at 30 eV penetrated into the inside of the graphite surface 
and were adsorbed between interlayer. The desorption of C2H2 on the clear graphite 
surface was observed in only the case incident at 5 eV. The animation of the MD sim- 
ulation and radial distribution function indicated that the graphenes were peeled off 
one by one at regular interval. In common to the incident energy the yielded molecules 
often had chain structures terminated by hydrogen atoms. The erosion yield increased 
compared with the case of no interlayer intermolecular force. 
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1 Introduction 

In the research into nuclear fusion, we deal with plasma surface interaction (PSI) prob- 
lem [1-7]. A portion of the plasma confined in an experimental device flows into a di- 
vertor wall, which is shielded by graphite or carbon fiber composite tiles. The incident 
hydrogen plasma which has weak incident energy erodes these carbon tiles in a pro- 
cess called chemical sputtering. The erosion produces hydrocarbon molecules, such as 
CH X and C2H X , which affect the plasma confinement. The PSI has been researched using 
molecular dynamics simulation (MD) [8-11]. 
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The authors investigated the PSI of graphite surfaces using the modified Brenner re- 
active empirical bond order (REBO) potential [12]. This MD simulation showed that if 
incident energy was 5 eV, many incident hydrogen atoms were absorbed by the graphite 
surface, while if the incident energy was 15 eV, most incident hydrogen atoms were re- 
flected. These behavior appear also in the case of deuterium and tritium incidence [13]. 
This absorption and reflection can be explained by the chemical reaction between a single 
hydrogen atom and a single graphene [14-16]. However, the number of absorbed hydro- 
gen atoms seems to be independent of this graphite erosion because although the hydro- 
gen atoms are absorbed by the first graphene on the surface only multiple graphenes 
are destroyed simultaneously. Because the incident hydrogen atoms pushed the graphite 
surface, the chemical bond between the first and second graphenes occurs. This was the 
trigger of the graphite erosion. 

However, the our previous work for PSI did not represent the interlayer intermolecu- 
lar interaction of the graphite because a suitable model was not known. For example, the 
existing model of the interlayer intermolecular interaction does not deal with "ABAB" 
stacking of the graphite structure. Before we look into the effect of the interlayer inter- 
molecular interaction, we should create new model of interlayer potential. 

We used MD simulation to investigate hydrogen atom sputtering on graphite surface. 
In S]2[ modified Brenner REBO potential and new model of interlayer intermolecular po- 
tential are denoted. We describe the simulation model in fj3] In ^4] we present simulation 
results. We discuss in Sj5[ This paper concludes with a 



2 Potential Models 

2.1 Modified Brenner REBO potential model 

We describe the model of Brenner reactive empirical bond order (REBO) potential [17] 
and our modification points. This potential model is created based on Morse potential 
[18], Abell potential [19] and Tersoff potential [20,21]. 
The potential function U is defined by 
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where r,y is the distance between the z-th and the y'-th atoms. The bond angle 9j ik is the 
angle between the line segment which starts at the z-th atom and ends at the j-th atom 
and the line segment which starts at the z-th atom and ends at the k-th atom, as follows: 
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where x^ = Xj — Xj is the relative vector of position coordinate from the j-th atom to the 
z-th atom, and ru is the distance between the z-th and the ;'-th atoms. The dihedral angle 
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B^jjf is the angle between the triangle formed by the j-th, the z-th and the k-th atoms and 
the triangle formed by the z-th, the j-th and the Z-th atoms. The cosine function of O^f 
given by 



is 
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The repulsive function VSi (ly - ) and the attractive function Vj^ (r^ ■) are defined by 



l + 9W) Auj ex P 



fk( r v)E B nlij]exP[-pn[ij] r . 



n=l 



(2.4) 
(2.5) 



The square bracket such as [ij] means that each function or each parameter depends only 
on the species of the z-th and the ;'-th atoms, for example Vq C , V^h and V^ H (= Vjj C ). The 
coefficients Q^, A^, oc^, B n ^ and jS„^ are given by TableU 

The cutoff function (r^ ) determines effective ranges of the covalent bond between 
the z'-th and the j-th atoms. Two atoms are bound with the covalent bond if the distance 
Tij is shorter than Dj^j 11 . Two atoms are not bound with the covalent bond if the distance r,y 
. The cutoff function f?.„ (r,y) connects the above two states smoothly 
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The constants D™.j n and C|^| x depend on the species of the two atoms (Table [4]). The 
cutoff function (r^-) distinguishes the presence of the covalent bond between the z'-th 
and the ;'-th atoms. 

The potentials VSi and in Eq. (|2.1|l generate two-body force, because both are the 
function of the only distance ru. The multi-body force is used instead of the effect of an 
electron orbital. In this model, by (t/},{# B },{0 DH }) in Eq. I|2.1[) gives multi-body force 
and is defined by 



1 r 
2 
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The first term i [• • •] generates three-body force except the effect of n electrons. The sec- 
ond term njy C in Eq. (|2.7|) represents the influence of radical energetics and n bond conju- 
gation [17]. The third term fcjj >H ({r},{0 DH }) in Eq. (|2.7[) derives four-body force in terms 
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of dihedral angles. These functions are composed of the production of cutoff functions 
fEj, (tij). Five- or more-body force are generated during chemical reaction. 

The function Z£ _7r ({r},{0 B }) in Eq. flZTJ is defined by 



qj-"({r},i<P})= 1+ E/fiflM^Ccos^W+^CAf, Nf)J 2 . (2.8) 

The function G, in Eq. 02.81 1 depends on the species of the z'-th atom. If cos$? )t >cos(109.47 o ) 
and the z'-th atom is carbon, G, is defined by 



Q(cos^)= [l-Q c (Mf)] G c (cos0? fc ) + Q c (M}hc(cos0? fc ). (2.9) 



If cosSL. < cos(109.47° ) and the z'-th atom is carbon, G, is defined by 



G ! (cos^) = G c (cos^). (2.10) 

And, if the z'-th atom is hydrogen, G, is defined by 

Q(cos^) = Gh(cos^). (2.11) 

Here Gq, 7c an d Gh are the sixth order polynomial spline functions. Though the spline 
function G, needs seven coefficients, the only six coefficients are written in Brenner's 
paper [17]. We determine the seven coefficients in table |5[ [6] and respectively. The 
function Q c and the coordination number M\ in Eq. \2.9\ are defined by 

1 (x<3.2), 
Q c (x) = { i[l+cos(27r(x-3.2))] (3.2<x<3.7), (2.12) 
(x>3.7), 



The constant A[«w in Eq. 02.81) is a weight to modulate a strength of three-body force, 
which depends on the species of the z'-th, the j-th and the k-th atoms. In comparison 
with Brenner's former potential [22], we set constants Ar^y as follows: 

Ahhh = 4.0, (2.14) 

Accc = A CC H = Achc = Ahcc 

= A H hc = Ahch = Achh = 0. (2.15) 

The function in Eq. (|2.8|) is required in the case that molecules forms solid struc- 
ture. The function Pm is the bicubic spline function whose coefficients depend on the 
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species of the z'-th and the j-th atoms (Table H). The parameters Njj and Njj are, respec- 
tively, the number of hydrogen atoms and the number of carbon atoms bound by the z'-th 
atom as follows: 
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The second term n^ c in Eq. (|2.7[) is defined by a tricubic spline function F^-j as 
where the variables are defined by 

NHE/[ft]('*), 



carbon carbon 
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with 



C N (x) : 



r 1 (z<2), 
i[l+cos(7r(x-2))] (2<x<3), 



(2.21) 







(x>3). 



The second and the third terms of the right hand of Eq. (|2.20|) are not squared. We 
note that they are squared in Brenner's original formulation [17]. By this modification, 
a numerical error becomes smaller than Brenner's formation. Table shows the revised 
coefficients for Fm . 



The third term b^ H ({r},{9 DH }) in Eq. (27) is defined by 



6f ({rUe^^T^N^) 



(2.22) 



where Tuj\ is a tricubic spline function and has the same variables as Fwi in Eq. (|2.18[) . The 
coefficients for is also revised due to the modified N™ n) (Table [TDl. In the present sim- 
ulation, the function Tr«i becomes Tec (2,2,5) f° r a perfect crystal graphene, and becomes 
Tec (2,2,3) or Tec (2,2,4) when a hydrogen atom is absorbed. 

The time step should be smaller than that of general CMD. To keep numerical error 
small, we set 5 x 10~ 18 s in the present simulation because the potential model has com- 
plex form by cutoff functions and spline function. 
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2.2 Interlayer intermolecular potential 

In the research for Interlayer intermolecular force, the binding energy has been well in- 
vestigated. However, experimental data are not nearly enough and ab-initio calculation 
cannot give us correct results yet [23]. Especially, information of the repulsion of the inter- 
layer is hardly reported. Therefore, now, we have no other choice to create the potential 
model artificially. 

We propose new interlayer intermolecular interaction potential model for graphene 
layers. First, simple intermolecular potential function between carbon atoms is defined 
by 



^(r) = A{V-(?- 1 ) 



(2.23) 



where r is the distance between two carbon atoms, n is the exponent of attraction, and 
A,dl,c are the parameters to determine binding energy. If n > a, the potential function has 
a local maximum of positive energy on r = c. Though we tried modelling the interlayer 
intermolecular potential, the challenge fell through. Figure llfA) shows a potential func- 
tion which consists of simple intermolecular potential of Eq. eq:simpleV. Such as this 
potential model, we hardly produce the difference of the potential minimum energy be- 
tween the three type of stacking of Fig. [TJ We consider that the difficulty comes from the 
use of only two body force. The attractive interaction is regard as two body interaction 
historically, such as Lennard-Jones potential, Morse potential and their combination Eq. 
eq:simpleV, because it is effective in the long range. However, the repulsive interaction 
is effective in the short range. The force in the short range are provided by chemical 
interaction. Therefore, the two body force is inadequate to approximate the repulsive 
force. 

The chemical interaction is generally represented by multi-body force in the MD sim- 
ulation. Here, we propose interlayer intermolecular potential using three body force. The 
product of the simple two body force V^L(r«) of Eq. 12.231 and special cutoff function Q/ 



gives us 
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The special cutoff function Cy depends on the angles between three atoms as 
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where r,y = r,-— fj, r^ = |r ; y| and cos0«£ = (rij-f^)/ The functions / a (cos0) are given 






O — Carbon atom and covalent bond of first layer 
O — Carbon atom and covalent bond of second layer 



Figure 1: The three type of stacking of the graphite 
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where c on = 0.25 and c Q ff = 0.35. The function (r) is equal to the cutoff function of the 
modified Brenner REBO potential: 
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where the parameters are denoted in the Table |2] Now, we set the parameters a. and c to 
keep the interlayer distance 3.35 A as follows: c=1.8 A and a=4.84. If ^=0.9961498,2.9884494 and 4.980749, 
the interlayer binding energy par atom becomes 20 meV, 60 meV and 100 meV, respec- 
tively (See Tabled]). Figure |2fB) show that the new interlayer intermolecular potential 
model provide the difference of the minimum potential energy between the three types 
of stacking of Fig. [TJ As a result, the structure of "ABAB" stacking Fig. [Ha) become the 
most stable state. 
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Table 1: The parameters of intermolecular potential. The parameters ^20,^60 ar| d ^.100 correspond to the 
coefficient A in VJlM and determine the binding energy of the interlayer par atom to 20 meV, 60 meV and 
100 meV, respectively 

n = 6 a = 4.84 c = 1.8 A 

A 20 = 0.9961498 eV A 60 = 2.9884494 eV A m = 4.980749 

Cnyi — 0.25 Cnti — 0.35 



Table 2: The constants for the cutoff function f^^( r ij)- They depend on the species of the ;-th and the j— th 
atoms. 
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Figure 2: The interlayer potential energy by using two body force only (A) and present new model (B). The 
symbols (a), (b) and (c) correspond to the three type of stacking of the graphite in Fig. [TJ 
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3 Simulation Method 

The graphite which consists of eight graphenes [24] and has "ABAB" stacking were set 
to the center of coordinates parallel to x-y plane. Each graphene consisted of 160 carbon 
atoms measuring 2.00 nm x 2.17 nm. The size of the simulation box in the x- and in- 
directions is equal to that of the graphenes with the periodic boundary condition. The 
inter-layer distance of the graphite was initially 3.35 A. The carbon atoms obeyed the 
Maxwell-Boltzmann distribution at 300 K, initially. During the simulation, two carbon 
atoms in the 8-th graphene were fixed to block the movement of whole of the graphite. 
One was the center atom of the 7-th graphene from the surface, and the other was located 
at the boundary of the 8-th graphene. The graphite surface was oriented to face the 
positive z-direction. 

Hundreds of hydrogen atoms were injected at regular time intervals of 0.1 ps parallel 
to the z-axis. The z-coordinate of the injection point was 60 A. The x- and y-coordinates 
of the injection point were set at random. The initial momentum vector (0, 0, po) was was 
defined by 

PQ^y/TmEu (3.1) 

where E\ is the incident energy, and m is the mass of the incident hydrogen atoms. 

We adopt NVE conditions, where the number of atoms, volume, and total energy 
are conserved, except for the addition of incident atoms and removal of outgoing atoms. 
The simulation time was developed using second order symplectic integration [25]. The 
chemical interaction was represented by the modified Brenner REBO potential. The inter- 
layer intermolecular interaction was represented by the new model in Sj2[ The interlayer 
binding energy is selected to 60 meV. To keep the computational error of total energy 
small, the time step was 5 x 10~ 18 s. 

4 Results 

We performed MD simulates for the three cases in which the incident energy of all hydro- 
gen atoms are set into 5 eV, 15 eV or 30 eV. In this section, simulation results are described 
with the story of PSI process. 

In the initial short time period, graphite surface were not broken. However, the differ- 
ence between the incident energy caused the difference of hydrogen atom adsorption on 
the graphite surface. Figure |3{a),IUa) and|5|a) show the snapshots of the MD simulation 
for PSI at t = 2.16 ps, at which more than 20 hydrogen atoms had done chemical inter- 
action with the graphite surface. From the figures, we noticed the amount of adsorbed 
hydrogen atoms and adsorption sites on the graphite surface. For incidence at 5 eV, a lot 
of hydrogen atoms were adsorbed by the graphite surface. The adsorption sites are the 
front of the first graphene, where graphenes are numbered from surface side. The posi- 
tive and negative side of each graphene in the direction of z are called front and backside, 
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Figure 3: The snapshot of the MD simulation for PSI in the case of the incident energy of 5 eV. Green and 
white spheres represent carbon and hydrogen atoms, respectively. 



respectively. For incidence at 15 eV, few adsorbed hydrogen atoms exist on the front of 
the first graphene. The animation of the MD simulation illustrated that hydrogen atoms 
except for the adsorbed one were reflected by the first graphene and went back to the 
positive direction of z. For incidence at 30 eV, a lot of hydrogen atoms were adsorbed be- 
tween the first and second graphene layers, that is, the backside of the first graphene or 
the front of the second graphene. A few hydrogen atom are adsorbed on the front of the 
third graphene. The animation of the MD simulation of incidence at 30 eV demonstrated 
the following dynamics. A lot of hydrogen atoms passed through a hexagonal opening of 
the first graphene, which is formed by six C-C bonds. After that, a half of them was ad- 
sorbed on the backside of the first graphene and the others flowed between the first and 
second graphene layers. When approaching the first or second graphene, the hydrogen 
atom was adsorbed. The hydrogen atoms which were adsorbed by the third graphene 
had penetrated the first and second graphene at a stretch. All of the hydrogen atoms 
which are not adsorbed by the graphite were reflected by the first graphene. Namely, the 
hydrogen atoms which penetrated into the inside of the graphite surface do not go out 
again. In addition, it is never seen that a hydrogen atom hits out a carbon atom and ejects 
it from a graphene, which is called physical sputtering. 

While the graphite surface maintained the graphene sheet structure, small hydrocar- 
bon molecules, such as CH X and CiHx, did not occur for the incident energy of 15 eV 
and 30 eV. However, for the incident energy of 5 eV, one C2H2 was generated keeping the 
first graphene flat (See Fig. |3{b)). After the atoms continued the above process, the first 
graphene was destroyed independent of the other graphenes. After that, from the second 
sheet, the graphenes were destroyed one by one with time. These destruction process is 



11 




Figure 5: The snapshot of the MD simulation for PSI in the case of the incident energy of 30 eV. Green and 
white spheres represent carbon and hydrogen atoms, respectively. 
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common to all cases of incident energy. 

To estimate the breakage of the graphite, radial distribution function was calculated. 
However, we cannot define the three dimensional volume in this simulation model be- 
cause there is no boundary in the z-direction. Two dimensional radial distribution func- 
tion g(r,t) of each graphene layer was defined. First, we defined n, (r,i) as the number of 
the carbon atoms which are located at a distance of less than r from the z'-th carbon atom 
at time t. The average n(r,t) is then given by 



where means summation in only one graphene. Consequently, the radial distribu- 

tion function is given by 



We calculated each graphene layer. We plotted the maximum values of the radial dis- 
tribution function g max (r,t) as a function of time (see Fig. [6). These maximum values 
always demonstrated the amount of the C-C bonds of length r = 1.42 Aand correspond 
to the number of sp 2 bonds. The decrease of g max (r,t) indicates the destruction of each 
graphene layer. As the incident energy increases, the speed of the decrease of g max (r,t) 
increases. It seems that the fast decrease of g max (r,t) of each graphene occurred at inter- 
vals. This was roughly regular interval. 

Because yielded molecules repeated chemical reaction, the species of the yielded molecules 
was not identified, In common, the yielded molecules had chain structures, for example 
C-C-C-C-H. Hydrogen atoms were often located at the edge of the chain molecules. To 
estimate erosion yield Y(t), we counted the number of the carbon atoms that moved to 
the region z > 24 A, where the first graphene is initially located on z = 11.7 A. Figure 
shows the erosion yield Y (t) as a function of time t. The erosion yield Y (t) increases with 
time t linearly. As the incident energy increases, the speed of the increase of Y(f) become 
faster and the yielded molecules started to be created earlier. As a result of the use of the 
interlayer intermolecular interaction, Y(t) increased. 

5 Discussion 

5.1 Initial short process 

In the initial short time period, the behavior of hydrogen atoms depended on incident en- 
ergy. This behavior can be explained by the research of the interaction between a single 
hydrogen atom and a single graphene because of the following three reasons. First, the 
graphite surface was maintained during the initial short time period. That is to say, after 
doing interaction with hydrogen atoms, the C-C bonds in the graphite are not broken 
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(a) Incident energy is 5 eV 



(b) Incident energy is 15 eV 
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Figure 6: The maximum value of the radial distribution function g max ( r ,t) for each graphene layers as a 
function of the time t. Only figure (d) is the result of the previous MD simulation which neglects the interlayer 
intermolecular force. 



(a) Hydrogen incidence at 5 eV 



(b) Hydrogen incidence at 15 eV 



(c) Hydrogen incidence at 30 eV 
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Figure 7: The erosion yields Y of carbon atoms as a function of time. 
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and the graphenes keep their structures. Second, because a hydrogen atom was reflected 
or adsorbed before the next incidence, only one hydrogen atom relates to one process. 
Third, the interlayer distance of the graphite are kept about 3.35 A. Because a hydrogen 
atom and a graphene start chemical (strong) interaction at a distance less than 1.6 A, in 
the graphite structure one hydrogen atom does not interact with two graphenes simulta- 
neously. We, therefore, discuss the behavior of hydrogen atoms in the initial short time 
period comparing with the research of the interaction between a single hydrogen atom 
and a single graphene. We had already researched the interaction between a single hy- 
drogen atom and a single graphene using MD simulation with modified Brenner REBO 
potential [14-16]. In that simulation, a hydrogen atom was injected into a graphene ver- 
tically. We classified interaction into three types, which is adsorption, reflection and pen- 
etration. Moreover backside adsorption was distinguished from front adsorption. Since 
the injection were repeated tens of thousands of times while changing incident position, 
we obtained the rates of three types interactions. The rates of the three type interactions 
depend on the incident energy as follows; If the incident energy is less than 1 eV, almost 
all of the interactions become the reflection due to pz-electron on the graphene surface. 
For the incident energy from 1 eV to 7 eV, the adsorption is dominant and has a peak rate 
at 5 eV. All of this adsorption is of the front of graphene surface. For the incident energy 
from 7 eV to 30 eV, the reflection is dominant and has a peak at 15 eV. As the incident 
energy increases from 15 eV, the rate of the penetration increases. For the incident energy 
of more than 30 eV, the penetration becomes dominant. Moreover, for the penetrate pro- 
cess, the hydrogen atom needs to expand a hexagonal opening of the graphene to pass 
through. If the incident energy is not sufficient to expand the hexagonal hole and to leave 
the graphene, the hydrogen atom are adsorbed on the front or backside of the graphene. 
As a result, the rate of the adsorption has a small peak around 25 eV. 

In the present MD simulation, the adsorption and reflection for the incident energy 
of 5 eV and 15 eV are derived from the incident energy dependence of types of the inter- 
action between a single hydrogen atom and a single graphene. In the incident energy of 
30 eV, we consider that the behavior of the hydrogen atoms is described by the combi- 
nation of several type of the interactions. The first interaction with the first graphene is 
similar to the interaction between a single hydrogen atom and a single graphene. When 
the hydrogen atom penetrates the first graphene, it reduces its kinetic energy. Therefore, 
the incident energy of next interaction with the second graphene shifts to the range for 
the reflection or adsorption. Since the kinetic energy was small for the penetration, the 
hydrogen atom stayed in the current interlayer region. However, because the loss of the 
kinetic energy due to the penetration depends on the incident point and timing, a few 
hydrogen atoms maintain and penetrate the second graphene. If the hydrogen atom is 
driven by higher energy, it seems to go into deeper region. 

In addition, the rate of interactions between a single hydrogen atom and a single 
graphene hardly depend on graphene temperature for the incident energy of more than 
1 eV. This fact supports the above consideration even if the hydrogen atom heat up the 
graphite surface. 
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5.2 Trigger of surface destruction 



We discuss the trigger of graphite surface destruction in this subsection. 

In the present MD simulation, small hydrocarbons, for example CH T and C2H T , were 
not created while the graphite surface maintain the graphene sheet structure. The de- 
struction of graphite surface seems to be melting or amorphization due to heat from the 
incident energy and adsorption energy. In the case of the incident energy of 5 eV, we 
observed only one chemical sputtering on the clean graphite surface, which is the C2H2 
desorption maintaining the first graphene layer flat (See Fig. 13b)). It is considered that 
the incident energy flux is too high for the chemical sputtering on the clean graphite 
surface to occur. However, when the interlayer interaction is used, no chemical sputter- 
ing was observed. Therefore, we have advanced one step toward real process to have 
introduced the interlayer intermolecular interaction. 

Here, we had investigated a graphene erosion process in hydrogen atom gas using 
MD simulation [26]. In this previous work, the four type of C-C bonds appears. If one 
carbon atom of a C-C bond is adsorbing a hydrogen atom, the C-C bond is called mono- 
overhang C-C bond. If both of the two carbon atoms of a C-C bond are adsorbing a 
hydrogen atom on the same side of the graphene, the C-C bond is called ortho-overhang 
C-C bond. If both of the two carbon atoms of a C-C bond are adsorbing a hydrogen atom 
on the opposite side of the graphene, the C-C bond is called para-overhang C-C bond. If 
no carbon atom of a C-C bond is adsorbing a hydrogen atom, the C-C bond is called flat 
C-C bond. This previous work demonstrated that the para-overhang C-C bond is the 
most breakable in the four type of C-C bonds. This fact does not also support the chemi- 
cal sputtering on the clean graphite surface because the para-overhang C-C bond hardly 
appears in the present MD simulation commonly to all case of incident energy. Of course, 
our simulation cannot represent thermal effects such as thermal desorption because sim- 
ulation time is so short. However, from point of view of chemistry also, a C-C bond is not 
broken by the thermal effects and a carbon of a graphene cannot adsorb more than two 
carbon atoms. Therefore, it is hard to create small hydrocarbon on the clean graphite sur- 
face. We consider that the chemical sputtering occurs on the graphite surface which has 
many defects, edge regions or amorphous regions rather than the clean graphite surface. 
The chemical sputtering on the carbon amorphous surface and the thermal desorption 
from a polymer were reported [8, 27] 

Here, we note that when the interlayer intermolecular interaction is not used, the 
trigger of the surface destruction was the covalent bonding between the first and sec- 
ond graphenes. This covalent bonding generated heat by using its binding energy and 
broke the graphene structure. In the present simulation, the covalent bonds between the 
graphene layers hardly appear. 
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5.3 Steady state of PSI 

The maximum values of the radial distribution function ? max 

(r,t) indicates the amount 

of sp 2 bonds. The decrease of g m ax(i"A) in Fig. [6] implies that the flat structure of the 
graphene is broken. As the interaction energy increases, the speed of the decrease of 
gmax(r,t) become faster and the start time when the hydrocarbon molecules are yielded 
became earlier. These facts, in relation to the previous subsection, indicate that the heat 
of the graphite surface were derived from the incident energy regardless of chemical in- 
teraction. The animation of the MD simulation showed that the graphenes were peeled 
off one by one from surface side. In addition, it imply this fact that gmax(r,t) of each 
graphene started to decrease at regular intervals in Fig. |6l In the previous work in which 
the interlayer intermolecular force is neglected, the hydrogen atoms pressed the graphite 
surface because of the absence of the repulsive force between the graphene layers. The 
first and second graphenes were bounded by covalent bonds. This covalent bonding 
generated heat by using its binding energy and broke the graphene structure simultane- 
ously. From comparison, it is considered that the interlayer intermolecular interaction 
played a important roll of the repulsion to resist the pressure due to hydrogen atom in- 
cidence. As a result, in the present MD simulation, the graphenes was not connected by 
covalent bonds and then they peeled off individually. This is a mechanism particular to 
PSI because in nano graphite material science, it is thought that the attraction part of the 
interlayer intermolecular force is important for making a molecular structure. 

The erosion yield Y(f ) increases with time t linearly. This linearity process is namely 
regard as steady state. The steady state is also adhered by the fact that the graphenes 
were peeled off at regular intervals. Of course, because the number of the graphite layers 
are finite, the steady state did not continue for a long time. 

The present MD simulation achieved the steady state without temperature control. 
Therefore, the present MD simulation perhaps differs from real PSI process. Though 
some temperature control methods control exists, the problem is not solved even if the 
methods is used. The temperature control methods usually brings about rapid cooling 
because we has to finish a cooling process in the MD simulation time, which is at most 
nano-seconds. If the thermostat for temperature control acts to the graphite surface di- 
rectly, the movement of the atoms are restricted. Consequently, the chemical interaction 
on the graphite surface become far from a real behavior. On the other hand, in the re- 
search of MD simulation, we often create or use potential models to achieve a realistic 
trajectory. Thereby, the speed of heat transport is also realistic, that is, so slow for the 
time scale of the MD simulation. As a result, if the thermostat is set to the region which 
is remote from the graphite surface, the heat cannot be transported to the thermostat. We 
have to create a new method of the temperature control in the near future. 

In the present work, because the covalent bonds between the graphene layers hardly 
occur, the heat of the surface is not transport to lower graphene layers. Therefore, the ero- 
sion yield Y(t) increased compared with the case in which the interlayer intermolecular 
force was not used. 
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6 Summary 

The new model of the interlayer intermolecular potential to represent "ABAB" staking 
of the graphite was proposed. We performed the MD simulation of the hydrogen atom 
sputtering on the graphite surface for the three cases of the incident energy of 5 eV, 15 eV 
or 30 eV. In the initial short time period, keeping the graphene structure flat, the hydro- 
gen atoms brought about the difference interaction process. The first graphene adsorbed 
a lot of the hydrogen atoms which were injected at 5 eV and reflected almost all of the 
hydrogen atoms which were injected at 15 eV. The hydrogen atoms which were incident 
at 30 eV penetrated under the first graphene and were adsorbed between the graphene 
interlayer. These process is similar to the interaction between a single hydrogen atom and 
a single graphene. The C2H2 desorption on the clear graphite surface was observed in 
only the case of the incident energy of 5 eV. However, the small hydrocarbon molecules 
except for the one C2H2 were not generated from the clear graphite surface. We discussed 
this fact by using the result of the research for the graphene erosion in hydrogen atom 
gas. The animation of the MD simulation and the maximum values of the radial distri- 
bution function g max (r,t) for each graphene indicated that the graphenes were peeled off 
one by one at regular interval. In common to there cases of incident energy, the yielded 
molecules often had chain structures which is terminated by the hydrogen atoms. The 
linear increase of the erosion yield Y(f) is regard as the steady state of the sputtering 
process. The erosion yield Y(t) increased compared with the case in which the interlayer 
intermolecular force was not used. 

Acknowledgments 

The authors acknowledge stimulating discussion with Dr. Arimichi Takayama. The nu- 
merical simulations were carried out using the Plasma Simulator at the National Insti- 
tute for Fusion Science. This work was supported in part by a Grand-in Aid for Ex- 
ploratory Research (C), 2007, No. 17540384 from the Ministry of Education, Culture, 
Sports, Science and Technology. This work was also supported by National Institutes 
of Natural Sciences undertaking for Forming Bases for Interdisciplinary and Interna- 
tional Research through Cooperation Across Fields of Study, and Collaborative Research 
Programs (No. NIFS07KDAT012, No. NIFS07KTAT029, No. NIFS07USNN002 and No. 
NIFS07KEIN0091). 

References 

[1] T. Nakano, H. Kubo, S. Higashijima, N. Asakura, H. Takenaga, T. Sugie, K. Itami, Nucl. 

Fusion 42 (2002) 689. 
[2] J. Roth, J. Nucl. Mater. 266-269 (1999) 51. 

[3] J. Roth, R. Preuss, W. Bohmeyer, S. Brezinsek, A. Cambe, E. Casarotto, R. Doerner, E. Gau- 
thier, G. Federici, S. Higashijima, J. Hogan, A. Kallenbach, A. Kirschner, H. Kubo, J. M. Layet, 



18 



T. Nakano, V. Philipps, A. Pospieszczyk, R. Pugno, R. Ruggieri, B. Schweer, G. Sergienko, 

M. Stamp, Nucl. Fusion 44 (2004) L21. 
[4] B. V. Mech, A. A. Haasz, J. W. Davis, J. Nucl. Mater. 241-243 (1997) 1147. 
[5] B. V. Mech, A. A. Haasz and J. W. Davis, J. Nucl. Mater. 255 (1998) 153. 

[6] A. Sagara, S. Masuzaki, T. Morisaki, S. Morita, H. Funaba, M. Goto, Y. Takamura, K. 
Nishimura, N. Noda, M. Shoji, H. Suzuki, A. Takayama, A. Komori, N. Ohyabu, O. Mo- 
tojima, K. Morita, K. Ohya, J. P. Sharpe, LHD experimental group, J. Nucl. Mater. 313-316 
(2003)1. 

[7] C. Garcia and J. Roth, J. Nucl. Mater. 196-198 (1992) 573. 

[8] E. Salonen, K. Nordlund, J. Tarus, T. Ahlgren, J. Keinonen, C. H. Wu, Phys. Rev. B 60 (1999) 
R14005. 

[9] E. Salonen, K. Nordlund, J. Keinonen, C. H. Wu, Phys. Rev. B 63 (2001) 195415. 
[10] D. A. Alman, D. N. Ruzic, J. Nucl. Mater. 313-316 (2003) 182. 

[11] J. Marian, L. A. Zepeda-Ruiz, G. H. Gilmer, E. M. Bringaand, T. Rognlien, Phys. Scr. T 124 
(2006)65. 

[12] H. Nakamura, A. Ito, Mol. Sim. 33 (2007) 121. 

[13] A.Ito and H. Nakamura, to be published in Thin Solid Films, (cond-mat/0709.2976). 

[14] A. Ito, H. Nakamura, J. Plasma Phys. 72 (2006) 805. 

[15] A. Ito, H. Nakamura, A. Takayama, submitted, l |cond-mat/0703377| ). 

[16] H. Nakamura, A. Takayama and A. Ito, to be published, (cond-mat/0705.3130). 

[17] D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni, S. B. Sinnott, J. Phys.: 

Condens. Matter 14 (2002) 783. 
[18] P. M. Morse, Phys. Rev. 34 (1929) 57. 
[19] G. C. Abell, Phys. Rev. B 31 (1985) 6184. 
[20] J. Tersoff, Phys. Rev. B 37 (1988) 6991. 
[21] J. Tersoff, Phys. Rev. B 39 (1989) 5566; 41 (1990) 3248(E). 
[22] D. W. Brenner, Phys. Rev. B 42 (1990) 9458; 46 (1992) 1948(E). 
[23] M. Hasegawa and K. Nishidate, Phys. Rev. B 70 (2004) 205431. 
[24] H. -P. Boehm, R. Setton, E. Syumpp, Pure & Appl. Chem. 66 (1994) 1893. 
[25] M. Suzuki, J. Math. Phys. 26 (1985) 601. 
[26] A. Ito and H. Nakamura, submitted. 

[27] M. Yamashiro, H. Yamada and S Yamaguchi, J. Appl. Phys. 101 (2007) 046108. 



19 



Table 3: The parameters for the repulsive function and the attractive function Vj^. They depend on the 
species of the i-th and the ;-th atoms. 



[if 



Parameter 



CC 

0.3134602960833 A 
10953.544162170 eV 

4.7465390606595 A -1 
12388.79197798 eV 
17.56740646509 eV 
30.71493208065 eV 
4.7204523127 A -1 
1.4332132499 A -1 
1.3826912506 A -1 



HH 



0.370471487045 A 
32.817355747 eV 

3.536298648 A -1 
29.632593 eV 
OeV 
OeV 
1.71589217 A -1 
OA" 1 
OA" 1 



CH or HC 
0.340775728 A 
149.94098723 eV 

4.10254983 A -1 
32.3551866587 eV 
OeV 
OeV 

1.43445805925 A -1 

OA" 1 
OA- 1 



Qlij] 



My] 



oc {i 
B 



m 

B 2 [ij] 

B m 

Hi) 



Table 4: The constants for the cutoff function (r;y). They depend on the species of the i-th and the ;-th 
atoms. 

m gg (A) p^)f x (A) 

CC 1.7 2.0 

CH 1.3 1.8 

HH 1.1 1.7 



Table 5: The parameters for the sixth order spline function Gcfcosf^). 



-1 -0.001 0.10400 

-1/2 0.05280 0.170 0.370 -5.232 

cos(109.47°) 0.09733 0.400 1.980 41.6140 

1 8.0 0.23622 -166.1360 



Table 6: The parameters for the sixth order spline function Yc(cos0?j.). 



cosel 


7c 


7c 


a 


7? 


cos(109.47°) 


0.09733 


0.400 


1.980 


-9.9563027 


1 


1.0 


0.78 


-11.3022275 
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Table 7: The parameters for the sixth order spline function Gn(cos8j jk ) . The parameters are determined under 



COS0B,=O. 



Parameter Value 

G H (0) 19.06510 
G^jo) 1.08822 

G£(0) -1.98677 
G^ 3) (0) 8.52604 

G$(p) -6.13815 

G^ 5) (0) -5.23587 

Gg } (0) 4.67318 



Table 8: Parameters for the bicubic spline function P[ij]{Njl ,Njj) . The parameters which are not denoted are 
zero. 





Value 


Pcc(l,l) 


0.003026697473481 


Pec (2,0) 


0.007860700254745 


^cc(3,0) 


0.016125364564267 


Pcc(l,2) 


0.003179530830731 


Pcc(2,l) 


0.006326248241119 


Pch(1/0) 


0.2093367328250380 


Pch(2,0) 


-0.064449615432525 


Pch(3,0) 


-0.303927546346162 


Pch(04) 


0.01 


Pch(0,2) 


-0.1220421462782555 


Pch(14) 


-0.1251234006287090 


Pch(2,1) 


-0.298905245783 


Pch(0,3) 


-0.307584705066 


Pch(1/2) 


-0.3005291724067579 
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Table 9: Parameters for the tricubic spline function F^y The parameters which are not denoted are 
zero. The function F^ satisfies the following rules: F^(N lr N 2 ,N 3 ) = F^ j ^(N 2 ,N 1 ,N 3 ), d Nl F^(N lr N 2 ,N 3 ) = 
att%(N 2 ,Ni,N 3 ), F [lj] (N 1 ,N 2 ,N 3 )=F [ij] (3,N 2 ,N 3 ) if Nj>3, and F [ij] (N 1 ,N 2 ,N 3 )=F [lj] (N 1 ,N 2 ,5) if N 3 >5, 
where d^. = d/dN;. 

Variables 



Function 


M 

N\ 


i\/2 




value 


7" 1 /at at at \ 

Fcc{Ni,N 2 ,N 3 ) 


1 


1 


-1 

1 


0.105000 


1 


1 


2 


—0.0041775 




1 


1 
i 


^ to R 

O LU J 


U.UIOUOJD 




9 
Z 


9 
Z 


1 
1 


H HQ/1 /I /IQCT 




Z 


z 


9 
Z 


U.U40OZOD1 




9 

Z 


Z 




U.U0U0OZ04 




2 


2 


4 


f) f)l 5441 1 7 

W . W 1 J 11 11/ 




z 


9 
z 




n n 




n 


1 


1 


04338699 







1 


2 


0.0099172158 







2 


1 


0.0493976637 




o 


2 


2 


-0.011942669 




o 


3 


1 to 5 


-0.119798935 




1 


2 


1 


0.0096495698 




1 


2 


2 


0.030 




1 


2 


3 


-0.0200 




1 


2 


4 to 5 


-0.030133632 




1 


3 


2 to 5 


-0.124836752 




2 


3 


lto5 


-0.044709383 


d Nl F C c(Ni r N 2 M) 


2 


1 


1 


-0.052500 


2 


1 


3 to 5 


-0.054376 




2 


3 


1 


0.0 




2 


3 


2 to 5 


0.062418 


9M 3 fcc(N 1/ N 2/ N 3 ) 


2 


2 


4 


-0.006618 


1 


1 


2 


-0.060543 




1 


2 


3 


-0.020044 


Fhh(Ni,N 2 ,N 3 ) 


1 


1 


1 


0.249831916 


Fch(Ni,N 2/ N 3 ) 





2 


3 to 5 


-0.009047787516128811 


1 


3 


lto5 


-0.213 




1 


2 


lto5 


-0.25 




1 


1 


lto5 


-0.5 


Table 10: Parameters for the tricubic spl 


ine function Tqq. 


The parameters which are not denoted are zero. The 


function Tqq satisfies the 


following rule: 


T CC (N V N 2 ,N 3 ) 


= T C c(Ni,N 2 ,5) if N 3 >5. 



Variables 



Function N\ N 2 N3 Value 



T C c(N lt N 2 M) 2 2 1 -0.070280085 

2 2 5 -0.00809675 



